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We study soliton dynamics in a system of two linearly coupled discrete nonlinear Schrodinger 
equations, which describe the dynamics of a two-component Bose gas, coupled by an electromagnetic 
field, and confined in a strong optical lattice. When the nonlinear coupling strengths are equal, 
we use a unitary transformation to remove the linear coupling terms, and show that the existing 
soliton solutions oscillate from one species to the other. When the nonlinear coupling strengths 
Q.^ j are different, the soliton dynamics is numerically investigated and the findings are compared to the 

, results of an effective two-mode model. The case of two linearly coupled Ablowitz-Ladik equations 

' is also investigated. 
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I. INTRODUCTION 



The study of discrete solitons is a central topic in nonlinear lattice dynamics [U, H, [1] . Of particular relevance is the 
robustness of soliton motion against perturbations and the possibility of controlling the soliton dynamics. Accordingly, 
many recent studies on this subject are focused on multi-component nonlinear lattices (see, e.g., [1, [1] and references 
therein). In this context, the inter-component coupling has a crucial role, as it may either destabilize the soliton 
propagation, or it can be used as a resource to move the soliton from one component to the other in a controllable 
manner. 

The interest in multi-component systems is considerably motivated by the experimental realization and control of 
mixtures of Bose-Einstcin condensates (BECs) composed either by different hyperfine states ji, T, 8] (including spinor 
^H^. BECs 0), or by different species [l^, |Tl|. Importantly, multi-component BECs can be trapped in strong optical 
lattices, as in the recent experiment reported in Ref. [12|. The optical lattice can be created by one or more pairs of 
counterpropagating laser beams, giving rise to a periodic potential in one, two, or three spatial dimensions [l3[. In 
J> ■ such a case, the dynamics of each (single-component) BEC is typically well described by a discrete dynamical model, 
lO namely the discrete nonlinear Schrodinger equation (DNLSE) - see also [13, E [Hi for reviews and references 
^""^ therein. Denoting by the BEC wavefunction on site j of the optical lattice, this model reads 
'nI" ■ 

o: 

^\ where K is the tunneling rate between neighboring sites, u is the nonlinear coefficient proportional to the s-wave 
• scattering length, and Vj accounts for an external potential that may be superimposed to the lattice. Notice that 
^ ' although Eq. ([1]) has been written for a one-dimensional (ID) setting, a generalization to higher-dimensions is 
. straightforward. The DNLSE is a typical example of a discrete dynamical system, whose properties have been 
' intensively studied [3, [H, [i3| • The interest for the DNLSE model arises from the fact that it has been successfully 
used to describe the dynamical properties of several systems (including, e.g., arrays of coupled optical waveguides 
[Tsj). In ID, the DNLSE is not integrable [13]; however, and as concerns the homogeneous case {Vj = 0), soliton- 
like wavepackets exist and can propagate for a long time as stable objects, as it can be shown, e.g., by variational 
approaches Il9l . l20ll . Furthermore, the dynamics of such traveling pulses has been investigated in detail in the literature 

[in, [M [13, [Mill- 
Based on the above discussion, here we consider the dynamics of a multicomponent BEC in an optical lattice 
described by a system of coupled DNLSEs, the different BEC components being different hyperfine levels [13|. Since 
the atoms of different components interact with each other, one has a nonZmear coupling between the different DNLSE; 
furthermore, one can induce a linear Rabi coupling by coupling the different hyperfine levels by an electromagnetic 
field [Hm. 

The prototypical system we consider is a two-component Bose gas in an external trapping potential: typically the 
condensates are different Zeeman levels of alkali atoms like ^^Rb. Experiments with a two-component ^^Rb condensate 
use atom states customarily denoted by |1) and |2); in particular, the states can be \F — 2,7np = 1) and |2, 2), like, 
e.g., in [13|, or |I,— I) and |2,I), like, e.g., in [SOI (see also the recent work [^). In general, the condensates |1) and 
|2) have different magnetic moments: then in a magnetic trap they can be subjected to different magnetic potentials, 
eventually centered at different positions and having the same frequencies (like in the setup described in [30| ) or 
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different frequencies [29j| ; in the latter work, the ratio of the frequencies of V2 and Vi is \/2- It is also possible to add 
a periodic potential acting on the two-component Bose gas [26| . 

The two Zeeman states |1) and |2) can be coupled by an electromagnetic field with frequency Woxt and strength 
characterized by the Rabi frequency Or. A discussion of (and references on) the experimental manipulation of 
multicomponent Bose gases are in [23, [2g] . 

In the case of a binary BEC mixture confined in a one-dimensional lattice the system of the two coupled DNLSEs 
takes the form 

Otpn'jj . \ 2 2 

= -K2 (V'(2)j-l +V'(2)j + l) +Ul2|l/'(l)jPV'(2)j +U22|V'(2)jPV'(2)j + V,-l/'(2)j + ■ (3) 

In Eqs. ©-([S]), 4'{a)j denotes the wavefunction of the condensate a [a — 1,2) on site j, Q.{t) is proportional to 
the Rabi frequency (i.e., the strength of the electromagnetic field), the nonlinear coefhcients Uap are proportional 
to the scattering lengths a^p for species a and species /3, accounting for intra- (a — 13) and inter- (a ^ (3) species 
interactions. In Eqs. ©-(HI) we also assumed that the external potential superimposed to the optical lattice is equal 
for both species. Furthermore, Kq. is the tunneling rate between nearest-neighbour sites for particles of the condensate 
a; for a two-component condensate, when the periodic potential is the same for both the hyperfine levels one has 
Ki — K2 , otherwise one has Ki ^ K2 jsj] . 

In this work we are interested in studying the dynamics of the system, which is prepared with a soliton initially 
present in a single component, after the linear coupling il is turned on. Our presentation is organized as follows. First, 
in Section II we show that for Ki = K2 it is possible to perform a unitary transformation to remove the linear coupling, 
analogous to the one valid for two linearly coupled continuous Gross-Pitaevskii equations js^ (see also earlier relevant 
work in the context of nonlinear optics [33|,[3J|); one can then use the results valid for a single-component DNLSE to 
determine the soliton dynamics. The case Ki ^ K2 is numerically studied, and our findings are compared in Section 
III with the results of a simplified two-mode model. In Section IV, we consider two linearly coupled Ablowitz-Ladik 
equations (which is an integrable variant of the DNLSE [sF]) showing results similar to those obtained for two coupled 
DNLSEs. Finally, in section V, we briefly summarize our findings and present some interesting directions for future 
work. 

II. EQUAL TUNNELING RATES 
A. General Lattice 

In this Section we study the dynamics of a soliton initially present in one of the components, when the linear 
coupling with the other component is turned on. Although we are interested in solitons in \D coupled DNLSEs, to 
show the generality of our approach we write the Eqs. ©-([S]) on a general lattice as 

= --ft^l +"12|V'(2)jPV'(l)j + K7"V'(l)j + f^(^)V'(2)j, (4) 

Q / 

= -K2 ^ ^(2)n + "l2|V'(l)jP'0(2)j + U22\'4'(2)j?i\2)j + ^jV'(2)j + (5) 

where the sum X]n(j) "-"'^ ^^'^ neighbors n of the site j. We will also consider equal interaction-strengths (un = 
U22 = U12 = u), as it is almost the case, e.g., for the two-component *^Rb BECs studied in [Hils^; in fact, in Ref. 
the ratios wn : uu : U22 = 0.97 : 1.00 : 1.03 were used, while in Ref. |2cj| un : uu : U22 = 1-00 : 1.00 : 0.97. 
When Ki = K2 Eqs. dl])-® can be rewritten in a more compact form as: 

ih^ = T^j + (V'/GV'j)^, + V,i:, + n{t)P^,, (6) 

where 
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and 

T^, = -KY,i^n. (8) 

Then, similarly to the case of the continuum counterpart of our model [3^. [sslls^l. we introduce the spinor field 0i 
through the relation 

V'^ = U{t)4>^, (9) 

where U{t) is given by 

^« = «.'HPXW1 = ( -^^^)). ,10, 

with I{t) = {l/h) Jl n{t')dt', and obtain, 

ih^ = + (0/G0,)^, + Vj^j, (11) 

which corresponds to Eq. ([6]) without the Rabi term proportional to ^{t)P. This transformation was originally 
introduced in Ref. [1^ (see also |3J])j where it was used to eliminate a constant linear coupling term by a change of 
polarization basis in a set of two coupled continuous NLS equations modeling optical pulse propagation in birefringent 
optical fibers, and later was used in the context of binary BECs in Ref. [33] ■ It is important to stress that the above 
transformation ((9l)- (|10[) eliminates the Rabi term only when G and P commute; in other words, the transformation 
is only possible in the so-called Manakov case (sH], where all interaction strengths are equal, i.e., uu = U22 = ui2- 
In the case of nonequal intra- and inter-species interaction strengths, i.e., un = U22 7^ U12, the transformation can 
still be applied but the final equation contains a more involved nonlinear interaction term. This physically relevant 
situation was recently considered in Ref. [37j . where the latter calculation was carried out for the continuous case. 

Let us now consider a soliton solution of Eq. ^ with the the linear coupling turned off (ri(t) = 0) with initial 
condition 

Ut-0)=[^'%^^^), (12) 

i.e., all particles are initially in the first component. Here, Vl'j(t) denotes the soliton solution of the single component 
DNLSE for the condensate 1, which we assume to be known. If at time t = tg the linear coupling is turned-on and 
the Rabi frequency has a general time dependence ri(t), then it holds 

V'.XM = 0.(^0) = (13) 

Due to the fact that Eq. pT|) . which describes the evolution of (p for < > to, is identical to the one describing the 
evolution of ip for t < to (i.e., without the Rabi coupling), it follows that for t > to one has: 



Using Eq. pH)) . one gets 



Equation (I15|) shows that the soliton can tunnel without losing its shape, and the effect of the transformation is just 
a change in the normalization. Hence, normalizing ^'j to 1 and denoting 

Na=J2\ V-W. I' (16) 
3 

the fraction of the number of particles in condensate a, one gets A^i = 1 and N2 — tor t < to, as well as 

Ni^cos^I{t), (17) 

and 7V2 = 1 — iVi for t > Iq. We notice that the result ^Tl} does not depend on the particular soliton solution chosen, 
nor on the momentum of the soliton. 



FIG. 1: Fraction of the number of particles iVi in the condensate 1 as a function of time for Q{t) — fio (solid line) and 
n(t) = f2o sin/3(t — io) (dotted line), where f3 = Q.o/h. In both cases, time is in units of Qq/H. 



B. ID Chain 



As an example of the previous results we considered the ID case, and in Fig. [l]we plot (for t > io) the number 
of particles A'^i in the first component, for a constant Rabi frequency Q{t) — Qq and for a sinusoidal Rabi frequency 
il{t) = flo sin /3{t — to) versus time. One can clearly see that in both cases iVi oscillates between the values one and 
zero, as expected. 

As a further application of the previous analysis we consider the ID homogeneous problem with Vj — 0. In this 
case, measuring time in units of h/2Ki and energies in units of 2Ki, the rescaled coupled DNLSEs read 

= - ^ i^wj+i + V'(ib--i) + A (iV'dbf + 1^(2), n v-d), + ^m{2)j, (18) 

- i^W + 4'(2)j-i) + A (I^Adbf + H(2)j\^) ^-(2), + ^muh^ (19) 

where A = u/2Ki, ui = U/2Ki and R = K2/K1 is assumed to be equal to 1. At t = all particles are supposed 
to be in the first component in a soliton-like wavepacket and lo — Q ior t < Iq. It is known that on a chain, the 
single-component DNLSE soliton-like solutions can propagate for a long time even if the equation is not integrable 
Q. Therefore, we consider, at t = 0, a gaussian wavepacket centered at ^(t = 0) = ^Oi with initial momentum pq 
and width ^{t = 0) = 70 = -v/cJo. Then, its temporal evolution can be analyzed by a variational approach, where we 
assume that the wavefunction has the following form 



(20) 



a 2 



where ^(t) and j{t) = -y/ a{t) are the center and the width of the density respectively, having conjugate momenta 
p{t) and d{t) (/C is a normalization factor). Writing the equations of motion for ^,p,a,6 one obtains p{t) = po smd 
^ = sinpoe^'', with rj = 1/ (2a) + aS'^ /8. Imposing the conditions 7 = and 5 = 0, it follows that for cospo < the 
value of A (when 70 ^ 1) for which a (variational) soliton solution exists is [14| : 

Aso^ = 2V^^-^^. (21) 
70 

The stability of the variational solutions for large times has been investigated numerically in [l^, . For the value 
Asoi given by Eq. (PT|) one has a{t) = ao, 6{t) — and ^ = sinpo • e~^/(^"°^ ~ sinpo- Notice that p = n corresponds 
to a gap soliton. 

If at t = to we turn on the Rabi switch with frequency uj{t), then we obtain 

m f cosX(t) . vl/)'(i) \ 

(,-*sinX(i).4j(<) j- ^22) 
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FIG. 2: Fraction of the number of particles A'^i in the condensate 1 as a function of time for p = n (solid line) and p — 37r/4 
(dashed line) obtained from the numerical solution of the coupled DNLSEs p8|l - (|19[ ) with the initial condition given by (|20p 
and A = Asoi- The two lines are indistinguishable and coincide with (|17|) . also plotted as a dotted line. We used the parameter 
values: A = 1, to = 10, cJ = 0.2, 70 = 50. 



The numerical solution of the coupled DNLSEs P^ - P^ . with the initial condition V'(i)j(^ = 0) = ^Jit = 0) given 
by (|20p and A — Agoi, turns out to be in excellent agreement with (|22|l for different values of p, as shown in Fig. [2] for 
uj constant. 



III. DIFFERENT TUNNELING RATES 



When the tunneling rates are different (i.e., R — K1/K2 ^ 1), the unitary transformation pUj) does not really 
simplify the problem. Actually, Eqs. (SI)-® can be written in the compact form ([6]) with the notation 



T-0J = TKy,K2-4^j = 

By performing the transformation (jlOp one finally gets 



-^2E„0) V'(2)n- 



(23) 



= '^kr.k,<t^j + {<l^3^Gcj>j)^j + V,<f>, - ayASit)C{t) ^ 



(24) 



where A = Ki - K2, S = sinJ, C = cosJ, Ki{t) = Ki ~ AS^{t), 7^2(0 = K2 + AS'^{t), and ay 



-i 

1 



It 



therefore follows that the effect of the transformation ^TU\i is to replace the general time dependence n{t) of the Rabi 
switch, by the last term of Eq. ([M]) . So, whenever the tunneling rates are different we study Eqs. (U)-® directly. 

In Fig. [3] we plot the time dependence of A^i for different values of R. It is observed that iVi does not reach the 
value zero, i.e., the soliton cannot be totally transferred to the other component. Moreover, it is clear that the larger 
the deviation of R from 1 is, the larger the minimum reached value of A'^i becomes, and therefore the more inefficient 
the transferring process becomes. 

These numerical results can be understood by means of a simple two-mode model, where for t > to the variational 
wavefunction is assumed to be of the form 



(25) 



Here, vE'j'^ is given by (j20[) and the variational parameters are now the center ^, the momentum p, the width ^/a and 
its momentum S, plus A^i, A'^2 (the number of particles in the two components) and the phases ipi, ip2- Note that a 
similar to Eq. (j25p variational ansatz has been used in the past to study soliton dynamics in two linearly coupled 
continuous NLS equations describing pulse propagation in dual-core optical fibers ^38i] . 
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FIG. 3; Fraction of the number of particles A'^i in the first component as a function of time for p — tt obtained from the numerical 
solution of the coupled DNLSEs for four different values of R, namely 0.5 (top-left), 0.9 (top-right), 1.5 (bottom-left), and 2.0 
(bottom-right); the rest of the parameters are the same as in Fig. (2] 



The Lagrangian C is given by 

dt 



^-^Ev^rl^r-^, (26) 



where the Hamiltonian Ti. is 



^ = - J H (^(lb-^(lb + l + C-C.) (^(2)jV^(2)j + l + C-C. 



- U! 



j 3 



A 



E (^(ib-'/'(2), + c-c-j + 2^(1 V'db- 1' + 1 V'(2), I' +2 1 Vdb- n ^(2), n ■ (27) 

j j 



Omitting the details, the equations of motion for p and ^ are found to be 

p = 0, (28) 

^ = (A^i + ii-iVa) sinp • e"", (29) 

(where t] — 1/ (2a) + a(5^/8), and those for a and 6 are 

(5 = (TVi + i?A^2) cosp • (-i - ,52) 2A/V^7^ (30) 

d = 2(iVi + i?iV2)a(5cosp • e"''. (31) 



One sees that for = 1, since A^i + A^2 = 1, Eqs. ((28)) - ((3T|) do not depend on the equations for A^i, N2, (fii, (p2- 
When i? ^ 1, the dynamics of the internal degrees of freedom is coupled with the phase- number dynamics. The 
equations for A^i and N2 are 



Ni = 2wV^/ViiV2sin(^, (32) 

N2 = -2ujyjNiN2Sin^p, (33) 

where the relative phase is given by 

ip^Lpi- ip2, (34) 

and evolves in time according to 

^ = [l - R) cosp ■ e-^^ - uj { - cos ^. (35) 
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FIG. 4: Fractional imbalance 2 as a function of time for p — tt obtained from the numerical solution (solid line) of the coupled 
DNLSEs for two different values of R, namely 1.5 (top), and 2.2 (bottom); the dotted line is the solution of Eqs. p7|l - (|38|l for 
the same values of R, namely 1.5 (top), and 2.2 (bottom). The values of the other parameters are the ones used in Fig[2] 



To provide an estimate for the behaviour of Ni_2 from the system (j28p - p5p . we introduce the approximation 
g-i/(2Q!) 2^ which is reasonable for broad sohtons. In this case, by introducing the fractional imbalance 

z = Ni-N2, (36) 

the following equations for z and ip are obtained 

i = —2lli\/i~z^ simp, (37) 

z 

ip = {1 — R) cosp + 2u! cos (p. (38) 

Vl — z'^ 

One readily realizes that Eqs. ([571) - ([55)) are the equations of a dimer [s^, but without the mass term (i.e., 
corresponding then to a non- interacting dimer) and with a detuning term. These equations describe the dynamics 
of a BEG in a double-well potential [4l|, \^ and analytical expressions for the quantities of interest can be found 
[39I [iol . [4^ . One of the effects of the detuning is to induce the deviations from the complete transfer of the particles 
from one mode to the other, which is what is numerically found. In Fig. [4l we compare the numerical results obtained 
from the coupled DNLSEs with the solutions of Eqs. ([37|) - ((38|) : it is clear that the agreement between the two is fairly 
good. 



IV. LINEARLY COUPLED ABLOWITZ-LADIK EQUATIONS 

The Ablowitz-Ladik equation is an integrable variant of the DNLSE [35| . which in dimensionless units can be 
expressed as 

Accordingly, a system of two linearly coupled Ablowitz-Ladik equations can be written as 

4 ( ) - ( : ^sr; ) + a o ^ + 1 + ( ) ^ m 

Notice that this version of the vector Ablowitz-Ladik equation for w = constitutes the integrable vector (Manakov- 
like) generalization of the one-component Ablowitz-Ladik model [l7| . However, it is different from other coupled 
models of this type with either nonlinear [i^ or linear [31 coupling. 

By defining (pi through the unitary transformation (jlOp , substituting in Eq. (|40p and multiplying by , one obtains 
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which is in fact the same as Eq. (|40)) but without the Rabi term. In obtaining Eq. (|4ip . one takes advantage of the 
fact that I i/'a)^ P + I V'(2b- P=l '/>(i)j P + I (t>{2)j 

Since the integrable Eq. (|4T|) (for which soliton solutions are analytically available [ij]) does not depend on w, 
the analysis presented in Section II is applicable and yields the same results for the dynamics of an (exact) soliton 
wavefunction that corresponds to all particles being initially in one of the condensates. 



V. CONCLUSIONS 



In this work we considered the soliton dynamics in a system of two discrete, linearly-coupled, nonlinear Schrodinger 
equations. We showed that there is a unitary transformation that can be applied to the system and has the effect of 
eliminating the time-dependent linear coupling terms in the case when the nonlinear coupling coefficients are equal. 
We showed that the solitonic solutions that describe the number of particles of a two-component Bose gas, can oscillate 
from one species to the other. In the case where the tunneling rates are different, although the transformation can 
still be made, it has no simplifying effect on the analysis. Hence, we resorted to a numerical study of the problem 
and found good agreement with the results obtained from an effective two-mode model. These results indicate that 
the efficiency of the transfer mechanism is reduced monotonically with the deviation from the equal tunneling rate 
limit (see also [37|). Finally, we showed that the same unitary transformation can also be applied in the analysis of 
a system of two linearly-coupled Ablowitz-Ladik equations, transforming the linearly coupled vector model into the 
well-known integrable vector Ablowitz-Ladik equation. 

There are many extensions that can be considered in connection with this work. The same considerations can be 
applied to higher dimensions, where it is known that while continuum nonlinear Schrodinger solitons are unstable to 
collapse, discrete ones may be stable for sufficiently weak tunneling p^. l45j. On the other hand, one can generalize 
the phenomena examined herein to the case of true spinor condensates, e.g. for spin-1 (or higher) bosonic systems 
of *^Rb or ^^Na which have been experimentally realized [i^ \^ and are under intense theoretical investigation 



48 , |49|, |50|, |51| . Notice that using a three- mode approach for analyzing the transfer would be a particularly relevant 
47| approach in that setting. 
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